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Abstract 

The all-female Caucasian rock lizard species Darevskia dahli and other parthenogenetic species of this genus reproduce 
normally via true parthenogenesis. Previously, the genetic diversity of this species was analyzed using allozymes, 
mitochondrial DNA, and DNA fingerprint markers. In the present study, variation at three microsatellite loci was studied in 
1 1 1 specimens of D. datiii from five populations from Armenia, and new information regarding clonal diversity and clone 
formation in D. dahli was obtained that suggests a multiple hybridization origin. All individuals but one were heterozygous 
at the loci studied. Based on specific allele combinations, 1 1 genotypes were identified among the individuals studied. 
Individuals with the same genotypes formed distinct clonal lineages: one major clone was represented by 72 individuals, an 
intermediate clone was represented by 21 individuals, and nine other clones were rare and represented by one or several 
individuals. A new approach based on the detection and comparison of genotype-specific markers formed by combinations 
of parental-specific markers was developed and used to identify at least three hybridization founder events that resulted in 
the initial formation of one major and two rare clones. All other clones, including the intermediate and seven rare clones, 
probably arose through postformation microsatellite mutations of the major clone. This approach can be used to identify 
hybridization founder events and to study clone formation in other unisexual taxa. 
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Introduction 

Naturally occurring unisexual reproduction is documented in 
less than 0.1% of all vertebrate species [1]. Among vertebrates, 
true parthenogenesis has been detected only in squamate reptiles, 
and especially among lizards [2] . Lacertidae was the first family in 
which this phenomenon was discovered [3]. All but one [1] of the 
parthenogenetic species originate from interspecific hybridization 
between bisexual species. They are normally diploid or triploid 
and commonly exhibit high levels of nuclear gene diversity 
because of fixed heterozygosity at codominant loci [4]. Each 
parthenogenetic species usually consists of several of clonal 
lineages, possibly caused by either mutation, multiple origins, or 
genetic recombination [5 — 7]. It is suggested that clonal diversity 
can also be correlated with the size of the area of origin [8], 
distinct ecological conditions of the habitat, and age of partheno- 
genetic species [9,10]. 

In 1989, Parker et al. summarized the aUozyme variation 
patterns observed among parthenogenetic species of several genera 
of lizards and proposed a hierarchical model for their origin, to 
predict whether they arose from single or multiple interspecific 
hybridization [10]. They concluded that a unisexual species with a 
single hybridization origin typically has littie variation and usually 
consists of a widespread common clone with a few rare clones. For 



species with multiple hybridization origins, the pattern of genetic 
variation is highly variable, witii random allele combinations. At 
the same time, the nature of the founder event at the origin of 
unisexual species can vary significantly. In some cases, high levels 
of genetic diversity in allozymes, combined with low variation in 
mitochondrial DNA (mtDNA), suggest multiple origins from a 
geographically restricted sample of females [11 — 13]. In other 
cases, in which there is low diversity for allozymes and mtDNA, 
the origins appear to have been strongly restricted, both 
numerically and geographically [14,15]. 

The genus Darevskia is of special significance, because these 
lizards have been the subject of extensive ecological and 
biogeographical studies [16], and because parthenogenesis has 
probably arisen several times within the group [17]. This genus 
includes 17 bisexual and seven truly parthenogenetic, diploid 
species of hybrid origin [16]. The allozyme patterns observed in 
the parthenogenetic Darevskia species are generally in accordance 
with Parker et al.'s model. Among these species, Darevskia dahli has 
a single origin [18]. 

D. dahli, which is one of the seven unisexual species, occurs in 
the southern Caucasus Mountains of Armenia and Georgia 
[16,17]. AUozyme and mtDNA data analysis identified the sexual 
parents of D. dahli as D. mixta and D. portschinskii [8,17]. At present. 
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D. mixta occurs only in Georgia, and D. portschinskii is distributed in 
southern Georgia, northern Armenia, and western Azerbaijan 
[19]. According to Darevsky et al. (1967), parental species are 
mainly allopatric [19]. In Armenia, D. dahli is topically distributed 
in forest regions at altitudes ranging from 700 to 1500 m [20]. It 
inhabits moderately dry bedrocks and gorge slopes in the forest 
zone, and also occurs among stony sites in mountain steppe 
habitats in formerly forested sites. Although D. dahli populations 
are rather numerous, they are usually isolated from one another. 
However, some of them can be sympatric with bisexual and/ or 
parthenogenetic lizards of the genus Darevskia. For instance, in the 
Dzoraget region (near Stepanavan), D. dahli is sympatric with the 
bisexual parental species D. portschinskii and with parthenogenetic 
D. armeniaca; near Sevkar village, D. dahli is sympatric with 
parthenogenetic lizards D. rostombekoii and D. armeniaca, and with 
bisexual hzards D. raddei. The interspecific hybrids identified in 
these zones are infertile triploids [21]. In the Papanino population, 
two distinct color varieties of Z). dahli are known [16], suggesting 
the possibility of two distinct clones in this taxon. They differ in 
belly coloration, which can be yellowish or bright yellow [20]. 
Some differences between these morphological forms have been 
observed regarding ecology and egg-number production. Bright- 
yellow lizards produce fewer eggs than do pale-yeUow lizards, and 
they tend to be less xerophilous [20]. Examination of allozyme 
variation in D. dahli revealed the existence of five clones; however, 
unlike what was observed for the color varieties, one clone was 
widespread and abundant, and four other clones were rare and 
geographically restricted [18]. Based on such apattern of allozyme 
clone distribution and the current allopatric occurrence of the 
parental bisexual species, these authors concluded that few 
bisexual individuals were involved in the formation of the D. dahli 
species, and that the initial successful clone was formed only once 
[IS]- 

In our DNA fingerprinting studies of Z). dahli, a rather high level 
of variation among individual fingerprint patterns was demon- 
strated when various microsatellite probes were used for blot 
hybridization [22,23]. However, the origin of this variation- 
whether it occurred as a result of point mutations or more complex 
genomic reorganization—remains obscure. New information re- 
garding genetic variation and clone formation in D. dahli can be 
obtained via molecular and genetic studies of individual genomic 
loci. 

The main objective of this study was to examine whether the D. 
dahli species originated from a single or multiple interspecific 
hybridization events. For this, we sequenced all allelic variants of 
three microsatellite loci found among 1 1 1 D. dahli spc( im(;ns from 
five Armenian populations, identified nucleotide variations among 
sequenced alleles, and estabhshed a distinct genotype for each 
individual. We hypothesized that the nucleotide variations found 
among alleles in the hybrid genomes of D. dahli were inherited 
from two parental species, and that their specific combinations 
formed by two alleles in genotypes may reflect independent 
hybridization events between bisexual individuals of parental 
ancestors. In this way, we detected 1 1 genotypes and identified at 
least three hybridization founder events in the D. dahli populations 
sampled. These data provide new knowledge regarding the 
evolutionary origin of D. dahli that suggest its multiple hybridiza- 
tion origin. Our approach based on the detection of genotype- 
specific markers formed by combinations of parental-specific 
nucleotide markers can be used to address similar problems in 
other unisexual taxa. 



Materials and Methods 

Parthenogenetic lizards were collected between 1997 and 2006 
from their natural habitats in Armenia. AU territories in which 
lizard samples were collected are governmental. There are no 

protected areas or private land on these territories. Therefore, no 
specific permissions were required for our sample collections in 
these locations. Our field studies did not involve endangered or 
protected species, and according to the Red Book of Armenia, the 
genus Darevskia was not counted as an endangered species in that 
period [20]. The study was approved by the Ethics Committee of 
Moscow State University (Permit Number: 24-01) and was carried 
out in strict accordance with their ethical principles and scientific 
standards. Blood samples were taken from the tail veins of Uzards 
under chloroform anesthesia, and all efforts were made to 
minimize suffering. DNA samples from lizard blood were isolated 
using a standard phenol-chloroform extraction method with 
proteinase K, and r("suspi'nded in TE buffer, pH 8.0. PGR 
amplification was performed using DNA samples from 1 1 1 Z). dahli 
individuals of five isolated populations in Armenia that were 
available to us: "Papanino" (» = 69), "Fioletovo" (n = 9), "Vaagni" 
(b= 17), "Dzoraget" (n = 6), and "Dendropark" [n= 10) (Figure 1). 
The loci Du215, Du281, and Du323 were amplified using the 
primer pairs described previously (Du215:5' — CAACTAGCAG- 
TAGCTCTCCAGA-3', 5' -CCAGACAGGCCCCAACTT- 
3'; Du281:5'-TTGCTAATCTGAATAACTG-3', 5'- 
TCCTGCTGAGAAAGACCA-3'; Du323:5'-AAGCA- 
GACTGTACAAAATCCCTA-3', 5' -ACTGATCTAAAGA- 
CAAGGTAAAAT-3') [24,25]. These primers were first devel- 
oped for the amplification of the respective loci in the 
parthenogenetic Uzard species D. unisexualis [24]. In the course of 
cross-species amplification, they were found to be suitable for 
homeological loci amplification in the congeneric species D. dahli. 
The Du215, Du281, and Du323 loci were chosen for genotyping 
analysis because each of them was polymorphic and contained a 
tetranucleotide microsatellite cluster that met the requirements of 
good resolution of the individual allelic PGR products in a native 
polyacrylamide gel electrophoresis, thus allowing the direct 
sequencing of the individual alleles. In general, the procedure 
for the isolation and sequencing of individual alleles was carried 
out as described previously [24,25]. PGR was performed on 50 ng 
of DNA in a total volume of 20 |a,L using a GenePak PGR Core kit 
(Isogene, Russia) and 1 |XM of each primer. The reaction 
conditions were: one cycle of 3 min at 94°C; 30 cycles of 1 min 
at 94°G; 40 s at the annealing temperature (58°C for Du215, 50°C 
for Du281, and 48°G for Du323); and 40 s at 72°C, followed by 
one cycle of 5 min at 72°C. PGR products (15 |a,L) were loaded 
onto an 8% native polyacrylamide gel (to separate allelic variants 
for each locus) and run for 12 h at 60 V. A 50 bp ladder 
(Fermentas) was used as a size marker. The amplification products 
were visualized by staining the DNA in the gel with ethidium 
bromid(;. The well-resolved individual PGR products, which 
corresponded to the two individual alleles of the locus, were 
excised from the gel, purified by ethanol reprecipitation, and 
sequenced directly in both directions using a chain termination 
reaction with an ABI PRISM BigDye Terminator ver. 3.1 on an 
Applied Biosystems 3730 DNA analyzer. This procedure was 
carried out twice for genotypes 5-11, because each of them was 
represented by only one individual. In all of these cases, allelic 
identity was checked and confirmed via the comparison of 
sequences obtained independently. In a separate experiment, the 
PGR product of the Du215 locus in genotype 1 1 was excised from 
the polyacrylamide gel, purified, and cloned into the pMos blue 
vector according to standard procedures (pMos Blue Blunt-ended 
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Cloning kit, RPN5110; Aniersham Biosciences). The clones were 
amplified in MOSBlue competent cells grown at 37°C, and 
sequenced as described previously [25] . 

The nucleotide sequences were compared using the MEGA ver. 
5.1 software [26]. All new data obtained from DNA sequencing 
were deposited in GenBank (FJ981 592 -FJ98 1598; FJ975076- 
FJ975081; FJ98 1 559 -FJ98 1561). 

The number of alleles, allelic richness, as a measure of allele 
counts adjusted for sample size, and expected heterozygosity, as a 
measure of gene diversity, were calculated per locus and per 
population using the Fstat software ver. 2.9.3.2 [27]. 

The TCS program implements a statistical parsimony approach 
for estimating genealogical relationships between D. dahli geno- 
types using the method of Templeton et al. (1992) [28]. This 
method provides estimates of gene genealogy at the population 
level that are more accurate than traditional phylogenetic methods 
[29]. It is more frequendy used for SNP comparison of maternally 
inherited mitochondrial genes, which are known for their 
neutrality and linear arrangement [30]. Usually, microsatellite 
repeats are also considered as evolutionarHy neutral DNA markers 
[31]. Stable reproduction of the maternal diploid number of 
chromosomes in successive generations caused by premeiotic 
chromosome duplication during oogenesis in parthenogenetic 
hzards allowed us to consider homeological alleles at one locus as 
different genes, i.e., presentation of the two alleles of a 
heterozygous genotype in a linear arrangement. We assumed that 
the processes of gain and loss of microsatellite repeat units have an 
equal probability. Genotypes 1—9 were used for the raw matrix, in 
which each GATA repeat was designated by the letter A (Table 
SI). The program collapses identical sequences into haplotypes, 
calculates the frequencies of the haplotypes in the sample, and 



connects them into a network. Haplotype networks were 
established using TCS software ver. 1.21, with gaps being 
considered as a second state. 

Results 

In total, 1 1 1 specimens of D. dahli from five populations of 
Armenia were analyzed using locus-specific PGR and DNA 
sequencing of PGR amplificants. It was shown that all D. dahli 
individuals but one were heterozygous at the three loci analyzed, 
and contained two alleles that differed from each other regarding 
the length and structure of microsatellite clusters, and regarding 
single nucleotide variations (SNVs) in fixed positions of the 
flanking regions (Table 1). The number of alleles varied from three 
to seven (depending on the locus). All allelic variants of every locus 
could be divided into distinct groups according to the fixed 
nucleotide variations in the microsatellite flanking regions. For 
example (Table 1), fixed nucleotide variations in alleles 1-4 of 
Du215 formed the SNV set TGC, whereas those in alleles 5 and 6 
of Du215 formed the SNV set ACC, and those in allele 7 of 
Du215 formed the SNV set AGT. 

To identify genotypic diversity in D. dahli populations, allelic 
combinations of the three loci in each of the 1 1 1 individuals were 
constructed. In total, 1 1 genotypes that differed in population 
frequencies and geographical distribution were revealed (Table 2). 
The individuals with identical genotypes formed distinct clonal 
lineages. One clone (genotype 1) was abundant (designated as 
major) and was represented by 72 individuals (64.9% of the total 
cohort) in three sampled populations; another (genotype 2) was less 
abundant (designated as intermediate) and was represented by 2 1 
individuals (18.9% of the total cohort) in four populations; nine 
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Figure 1 . Map of Armenia with the distribution of localities from which parthenogenetic Darevskia tya/?// populations were collected. 

Sampling localities are indicated by different colours. Papanino (+40° 42' 27.76", +44° 45' 43.89"), Phioletovo (+40° 44' 29.53", +44° 39' 31.25"), 
Vaagni (+40° 52' 19.16", +44° 31' 28.60"), Dzoraget (+40° 54' 15.02", +44° 40' 37.78"), Dendropark (+40° 59' 30.55", +44° 21' 50.79"). 
doi:1 0.1 371/journal.pone.0091 674.g001 
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Table 1. Allelic variations of microsatellite containing loci Du215, Du281, and Du323 in parthenogenetic D. dahli species. 



Allelic 
vsridnt 


Size 

(bp) Structure of mlcrosatelilte cluster 


Fixed nucleotide 
variations^ 




EMBL/Gene 
Bank ac> No> 


Source 


Du215(dahli)l 


252 


5' 


(GAT){GACA){GATA)io(GACA)7{GATA){GCAA) 3' 


T (-58), G (-38), C (- 


-19) 


FJ981592 


this study 


Du215(dahli)2 


248 


5' 


(GAT){GACA)(GATA)9 (GACA)7(GATA){GCAA) 3' 


T (-58), G (-38), C (- 


-19) 


FJ981593 




Du215(dahli)3 


244 


5' 


(GAT){GACA){GATA)s {GACA)7(GATA){GCAA) 3' 


T (-58), G (-38), C (- 


-19) 


FJ981594 




Du215(dahli)4 


232 


5' 


(GAT)(GACA)(GATA)7 (GACA)5(GATA){GCAA) 3' 


T (-58), G (-38), C (- 


-19) 


FJ981595 




Du215(dahli)5 


232 


5' 


(GAT){GATA)ii{GCAA)4 3' 


A (-58), C (-38), C (- 


-19) 


FJ981596 




Du215(dahli)6 


228 


5' 


(GAT)(GATA)io(GCAA)4 3' 


A (-58),C (-38), C (- 


19) 


FJ981597 




Du215(dahli)7 


192 


5' 


fGATVGATAU 3' 


A (-58), C (—38), T (- 


-19) 


FJ981598 




Du281 (dahli)l 


225 


5' 


(GGTA){GATA)9(GAT)(GATA)(GGTA)2(GAT)(GATA)4 3' 


T (+15) 




FJ975076 


Davoyan et al. [42] 


Du281(dahli)2 


221 


5' 


(GGTA)(GATA)8(GAT)(GATA)(GGTA)2(GAT)(GATA)4 3' 


T (+15) 




FJ975077 


this study 


Du281(dahli)3 


199 


5' 


(GATA),2 3' 


C (+15) 




FJ975078 


this study 


Du281(dahli)4 


195 


5' 


(GATA),, 3' 


C(+15) 




FJ975079 


Davoyan et al. [42] 


Du281(dahli)5 


191 


5' 


(GATA),o 3' 


C (+15) 




FJ975080 




Du281(dahli)6 


183 


5' 


(GATA)8 3' 


C(+15) 




FJ975081 


this study 


Du323(dahli)1 


215 


5' 


(AOs. . .(GATA)i 1 (GAT)(GATA)(GATATAT){GA)4 3' 


C (-23), T (+39) 




FJ981559 


Vergun et al. [43] 


Du323(dahli)2 


211 


5' 


(AC)6...{GATA)io(GAT)(GATA)(GATATAT){GA)4 3' 


C (-23), T (+39) 




FJ981560 




Du323(dahli)3 


184 


5' 


(AC)5...{GATA)(GGT){GATA)3(GAT){GATATAT)(GA)4 3' 


A (-23), C (+39) 




FJ981561 





^Distances before ( — ) and after (+) microsatellite cluster are given in bp. 
doi:l 0.1 371 /journal.pone.0091 674.t001 



clones were rare, geographically restricted, and were represented 
by one (genotypes 5-11; each in 0.9% of the total cohort) or six 
and five (genotypes 3 and 4, 5.4% and 4.5% of the total cohort, 
respectively) individuals. 

One can see that the genotypic diversity among the five lizard 
populations varied from 0% to 44.4% (average, 9.9%) (Table 2). 
The highest levels of genotypic diversity were observed in the 
Vaagni and Dendropark populations. We calculated the number 
of population genetic indices for k loci using the Fstat software 
(Table S2). The genetic variations of genotypes 1—9 were used. 
The estimates of expected heterozygosity varied from 0.5 1 to 0.65 
(average, 0.56). The number of alleles varied from 2 to four with 
2.7 in average. The values of allelic richness varied from 2.000 to 
3.453 with 2.570 in average. The highest aUehc scores, the number 
of alleles as well as allelic richness (Rs) were observed in two 
populations for the aU three loci Dendropark and Vaagni: 3.453 
and 2.938, respectively (Du215), 2.600 and 2.800, respectively 
(Du281), 2.853 and 2.648, respectively (Du323). 

To clarify whether the D. dahli species originated from single or 
multiple interspecies hybridization events, genotype-specific nu- 
cleotide markers formed by combinations of parental-specific 
SNVs of three loci were established and compared. Figure 2 
schematically shows the structural composition of 1 1 genotypes 
that included microsatellite clusters and parental-specific SNVs in 
fixed positions in the flanking regions (Table 1). It was clear that 
genotypes 1-9 had identical combinations of parental-specific 
SNVs for aU three loci: TGC/ACC (Du215), T/C (Du281), and 
CT/AC (Du323). This means that they have identical genotype- 
specific markers and, therefore, may have a common origin as the 
result of a single interspecies hybridization event. However, the 
possibility of independent crossings of the parental individuals with 
frequent respective allelic variants for some of these genotypes 
cannot be ruled out. Genotypes 1 0 and 1 1 differ from each other 
and from genotypes 1-9 by parental-specific SNV combinations 
for locus Du215: TGC/ACT and TGC/TGC, respectively. Thus, 
we suggest that genotypes 10 and 11 have independent origins. 



Consequendy, at least three independent interspecific hybridiza- 
tion events took place in the genesis of the 1 1 genotypes. Several 
additional lines of evidence suggest that the microsatellite clones 
observed in D. dahli result from at least three independent 
interspecies hybridization events, rather than from a single origin. 
In addition to having the same combination of parental-specific 
SNVs, genotypes 1-9 also exhibited a similar organization of 
allelic microsatellite clusters. Different situations occurred in 
individuals with genotypes 1 0 and 1 1 . The microsatellite cluster 
of the ACT allele of genotype 1 0, unlike what was observed for all 
other alleles, had no (GCAA)4 nucleotide group. The individual 
with genotype 1 1, unlike all other individuals, was homozygous at 
the Du215 locus. 

An unresolved question remains regarding the relationships 
between genotypes 1-9. These genotypes differed from each other 
only by microsatellite sequences. Figure 3 schematically shows tiie 
resulting SP network, which reflects the relationships between 
genotypes 1-9 (genotypes 10 and 11 are plotted separately). The 
areas of the circles are correlated with the abundance of the 
genotype. The areas of the circle segments are correlated with the 
population distribution for each genotype. The black circles show, 
unsampled but computer-predicted genotypes that participated in 
the formation of genotypes 5, 6, and 8. 

According to Figure 3, the most distant genotypes were unique 
genotypes 6 and 5 from Papanino and Vaagni, respectively. It 
should be noted that the maximum variability of genotypes was 
observed in the latter population (genotypes 1, 2, 4, 5, and 8—11). 
Three of them (5, 8, and 9) were specific for the Vaagni 
population. Two genotypes (3 and 7) were specific for the 
Dendropark population among the three genotypes discovered. 
Only genotype 6 of the four genotypes detected was specific in the 
Papanino population, and no specific genotypes were found in the 
Fioletovo and Dzoraget populations. Genotypes 4, 1, and 2 were 
identified in two, three, and four populations, respectively. The 
analysis of the spatial-frequency distributions of the three more 
widely distributed genotypes and of a few population-specific 
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Du215 



Du281 



Du323 



TQC 
ACC 
TOC 

Arc 

TGC 
ACC 
TGC 
ACC 



TGC 
ACC 



TGC 
ACC 



TGC 
ACC 



TGC 
ACC 



TGC 
ACC 



10 

11 

I TGC 



(QAT) (GACA) (GATA)b(QACA)7(0ATA) (OCAA) 

■ (OAT) (aATA)iD (aCAA]4l 
(GAT) (GACA) (GATA)B(GACA)7(GATA) (GCAA) 

■ (OAT) (OATA)io (QCAA]<I 
(GAT) (GACA) (GATA)9(GACA)7(GATA) (GCAA) 3 

■ (GAT) (OATA)to (GCAA)*! 
(GAT) (GACA) (0ATA)b(GACA)7(GATA) (GCAA) ^ 

■ (GAT) (OATA)ii (GCAA)*! 
(OAT) (GACA) (OATA)8(GACA)7(GATA) (GCAA) 

■ (GAT) (OATA)ii (GCAA)*^ 
(GAT) (GACA) (GATA)b(GACA)7(GATA) (GCAA) 

■ (OAT) (OATA]io (GCAA)4l 
(OAT) (OACA) (OATA)io(QACA)7(QATA) (OCAA) 

■ (OAT) (OATA)iD (GCAA)4l 
(GAT) (GACA) (QATA)9(GACA)7(GATA) (GCAA) 

■ (GAT) (aATA)iD (GCAA)*! 
(OAT) (OACA) (0ATA)B(GACA)J(GATA) (GCAA) 

■ (GAT) (GATA)iD (GCAA)«! 
(OAT) (OACA) (0ATA)7(GACA)s(0ATA) (GCAA) 
(OAT) (OATA)s I 
(OAT) (GACA) (0ATA)8(QACA)7(GATA) (GCAA) 
(OAT) (GACA) (0ATA)8(GACA)7(QATA) (GCAA) I 



■ (GOTA) (QATA)9(QAT) (OATA) (GGTA);(OAT) (QATA)<| 

(GATA)ii 

■ (GGTA) (GATA)9(GAT) (GATA) (GGrA)!(GAT) (GATA)*! 

■ (OATA) u 

■ (GOTA) (GATA)9(0AT) (GATA) (OGTA);(GAT) (GATA)«! 

■ (GATA) 11 

■ (OCITA) (GATA)9(GAT) (GATA) (GGTA)z(OAT) (GATA)<I 

■ (OATA) 11 

■ (GGTA) (GATA)b(GAT) (GATA) (GGrA)!(GAT) (GATA)*I 

(GATA) 11 

■ (GOTA) (GATA)9(GAT) (GATA) (GOTA)z(GAT) (GATA)*1 

(GATA)e 

i(GOTA) (GATA)9(GAT) (GATA) (GOTA)!(GAT) (GATA)*I 
(GATA) 11 

■ (GOTA) (GATA)9(GAT) (GATA) (GOTA)!(GAT) (GATA)*! 

(OATA) 10 

■ (OOTA) (0ATA)9(0AT) (OATA) (GOTA)2(0AT) (OATA)*l 

(OATA) 12 

■ (GOTA) (GATA)9(GAT) (GATA) (GOTA)!(GAT) (GATA)*! 

(OATA) 8 

■ (OOTA) (0ATA)9(0AT) (OATA) (0OTA)2(0AT) (OATA)*l 

■ (OATA) 11 



I (AC)i...(QATA)n(QAT)((UlTA] (6ATATAT) (OA) 
(AC)s. . . (GATA) (GOT) (OATA) J (GAT) (GATATAT) (GA) 
(AC)4. . . (GATA) lO(GAT) (GATA) (GATATAT) (GA) 
(AC)s. . . (GATA) (GGI) (GATA)3(GAT) (GATATAT) (GA) 
I (AC)i...(aATA)ii(OAI)(aATA) (GATATAT) (GA) 
I (AC)s. . . (GATA) (GOT) (GATA)] (GAT) (GATATAT) (GA) 
I (AC)»...(GATA)io(GAT)(GATA) (GATATAT) (GA) 
I (AC)s. . . (GATA) (GOT) (OATA))(OAT) (GATATAT) (OA) 
I (AC)i...(QATA)io(GAT)(GATA) (GATATAT) (OA) 
I (AC)s. . . (GATA) (001) (OATA)3(0AT) (OATATAT) (OA) 
I (AC)6...(0ATA)ii(GAT)(GATA) (GATATAT) (OA) 
I (AC)s. . . (OATA) (GOT) (QATA)3(GAT) (GATATAT) (GA) 

■ (AC)(. . . (QAIA)ii(OAT) (OATA) (OATATAT) (OA) 

■ (AC)s. . . (OATA) (001) (OATA)i(OAT) (OATATAT) (OA) 
I (AC)e...(GATA)ii(GAT)(QATA) (OATATAT) (OA) 
I (AC)J. . . (OATA) (OOT) (OATA)}(QAT) (OATATAT) (OA) 

■ (AC)». . . (GATA) lO(GAT) (GATA) (GATATAT) (GA) 

I (AC)s. . . (GATA) (GOT) (GATA)}(GAT) (GATATAT) (OA) 
I (AC)&. . . (GATA) iD(GAT) (GATA) (GATATAT) (OA) 
I (AC)s...(GATA) (GGr)(GATA))(GAT)(GATATAT)(GA) 
I (AC)6. . . (OATA) lo(OAT) (GATA) (GATATAT) (OA) 
I (AC)s...(GAIA) (GOT) (GATA) 3 (GAT) (GATATAT) (GA) 



Figure 2. Schematic representation of 11 genotypes formed by allelic combinations of microsateilite loci Du21 5, Du281, and Du323 

in 111 D. (ya/)// individuals. Parental-specific SNV markers are shown in square. Variable microsateilite clusters are shown in each of two alleles. 
doi:10.1371/journal.pone.0091674.g002 



genotypes revealed a dependence between the frequencies of the 
clones and the geographical disjunction between the Vaagni 
population and two groups of populations, Papanino/Fioletovo 
and Dzoraget/Den(iropark. Namely, genotypes 1, 2, 4, and 6 were 
(discovered in the Papanino and Fioletovo populations, which are 
located at a distance of 9.5 km from each other. Genotypes 2, 3, 
and 7 were distributed in the Dzoraget and Dendropark 
populations, which are located at a distance of 28 km from each 
other. It follows that genotype 2 was shared by the two population 
groups. First, in the Vaagni population, which is located closer to 
the Dzoraget/Dendropark group (13— 19 km) than it is to the 
Papanino/Fioletovo populations (18 — 27 km), the common geno- 
type 2 was present, as were genotypes 1 and 4, which were also 
discovered in the Papanino/Fioletovo populations. Second, 
genotype 8 was discovered only in the Vaagni population and 
can be regarded as a derivative from genotype 3, which was 
specific to the Dzoraget/Dendropark group. Thus, the Vaagni 
population might be the source of the clonal diversity observed 
among the five populations studied. 

Discussion 

The main finding of this study was the presence of at least three 
interspecific hybridization founder events in D. dahli, which 
resulted in the formation of three initial clones and eight clones 
that probably arose through postformation microsateilite muta- 
tions of the initial clones. Previously, an aUozyme analysis 
conducted at 35 gene loci in three Armenian and three Georgian 
sampled populations of D. dahli revealed one major widespread 
clone and four rare clones detected only in restricted populations 
[18]. The authors concluded that the observed allozyme pattern in 
D. dahli generally followed Parker et al. 's model [10], which is 
consistent with its single hybridization origin. Nevertheless, they 
noted that possible multiple origins for D. dahli should be 



investigated further using more effective genomic markers [18]. 
Among parthenogenetic Darevskia, D. rostombekovi was also suggest- 
ed to be of single origin [32] . In D. armeniaca, four allozyme clones 
were detected [33], but one rare clone made up the majority of 
two populations, and another rare clone had different allelic 
compositions at two loci. The authors assumed that both mutation 
and multiple origin could be alternatively considered as a possible 
explanation of the origin of the clonal variation in D. armeniaca. 
However, the allozyme data failed to distinguish multiple origins 
and mutation events. Fu et al. (2000) noted that many studies that 
were performed at the species level examined various population 
characteristics, but few did this at the genotype level [33]. 

In our study, a new approach for the detection of D. dahli 
genotypes formed by independent hybridization between bisexual 
parental individuals was developed. This approach included 
microsateUite genotyping of the sampled individuals, and the 
identification of genotype-specific nucleotide markers formed by 
combinations of parental-specific nucleotide markers (SNVs). The 
examination of the allelic combinations of three microsateUite loci 
in 1 1 1 D. dahli individuals revealed 1 1 clones, two abundant (major 
and intermediate) and nine rare and geographically restricted. 
This pattern of microsateUite-clone distribution raises questions 
about the scenario of the D. dahli species formation, and about the 
relationship between major, intermediate, and rare clones. The 
comparison of genotype-specific markers showed that all individ- 
uals with genotypes 1-9 exhibited the same combination of 
parental-specific SNVs and, therefore, could have a common 
origin. Two individuals with genotypes 10 and 11, respectively, 
possessed a unique combination of parental-specific SNVs, and 
thus were likely to have independent origins. Our data do not rule 
out further hybridization events. First, the variation observed 
among clones 1 — 9 could actually be a polymorphism that was 
inherited from the parental populations, and there were differ- 
ences between individuals with similar genotypes; that is one 
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Figure 3. Schematic representation of the SP networi< that reflects relationship between genotypes 1-9 in D. dahli. Complete 
sequences of D. dahli genotypes were analyzed using TCS software version 1 .21. Genotypes 10 and 1 1 are plotted separately. Population distribution 
of the genotypes is shown by different colours. The black circles show, unsampled, but computer-predicted genotypes. 
doi:10.1371/journal.pone.0091674.g003 



potential source of underestimation of the frequency of hybrid- 
ization. Second, in tlie case of low variation in parental 
populations, some hybridization events could have led to identical 
unisexual genotypes with fewer hybridization events. Nevertheless, 
some of the clones 1—9 may indeed have originated via 
postformation microsatelUte mutations of the initial clones. 
Consideration of the SP network, which reflects the relationships 
between genotypes 1-9, does not allow the disclosure of the initial 
clone. According to Parker et al.'s model [10], we predicted that 
the major widespread clone with genotype 1 could be the initial 
clone, whether all others probably arose through postformation 



microsatellite mutations of the major clone. Figure 3 shows that 
from one to four mutational events led to the formation of new 
microsatellite clones. Previously, we characterized the highly 
unstable (GATA)„-containing locus in parthenogenetic D. unisex- 
ualis. We identified and sequenced mutant alleles in some offspring 
of the first generation and showed that de novo mutations occurred 
via deletion or insertion of a single microsatellite repeat [25]. 

Multiple hybridization events leading to multiple aUozyme 
clones have been described in Cnemidophorus and Heteronotia [34]. 
The patterns of allozyme variability among species of multiple 
origins ranged from a large variation {Heteronotia hinoei; [6]) to 



PLOS ONE I www.plosone.org 



7 



March 2014 | Volume 9 | Issue 3 | e91674 



Clonal Diversity and Clone Formation in D. dahli 



medium [Cnemidophoms tesselatus; [35]), and slight {Darevskia 
parthenogens; [18,32]) variations. Thus, the patterns of aHozyme 
variation in D. dahli [18] do not contradict the multiple 

hybridization origin demonstrated in our study, and the higher 
microsatellite clonal diversity reflects the higher efiectiveness of the 
markers used. 

Nine relatively rare clones, each represented by one or several 
individuals, were found in four sampled populations of Z). dahli. Of 
the rare clones, six, including those with genotypes 1 0 and 1 1 , 
occurred in the Vaagni population. This population, together with 
the Dendropark population, has an enhanced level of clonal 
diversity (Table 2, Figure 3). Interpopulation differences formed by 
the distribution of various genotypes can be explained by the 
distinct ecological conditions of the habitat of D. dahli populations. 
However, differences in suitable landscapes, rock dimensions, and 
the resource utilization cur\'es are secondary for D. dahli with 
respect to the differences regarding preferred climates [36]. This 
species has a narrow climatic niche and its distribution depends 
mainly on both temperature and humidity. Thus, the varied 
topography, vegetation, and climate can influence some disjunct 
distribution and interpopulation differences of Z). dahli populations. 
Some other factors, which are connected with the frequency of 
interspecific mating at some locEilities and hybrid fitness ability, 
could also be significant in the case of the sympatric existence of 
the parental bisexual species in the past. At present, only one of the 
studied D. dahli populations (Dzoraget) overlaps with the bisexual 
parental species D. portschinskii, but the upper altitudinal limit of the 
unisexual species is higher [19]. The interaction between D. dahli 
and D. portschinskii is restricted to competition for food or cover 
and, as noted above, the interspecific hybrids observed were 
infertile triploids [21]. 

The greater clonal diversity of the Vaagni and Dendropark 
populations could be also connected with their geographical 
location in the seismically active region of Armenia, where a 
casting of chemical mutagens has sometimes occurred [37]. Under 
these conditions, the accumulation of microsatellite mutations 
might be increased. In previous studies, a high level of 
chromosomal instability was observed in some rodent and lizard 
species from this region, as well as chromosomal speciation events 
[38-41]. 

The lack of statistical approaches for the evaluation of 
population diversity in parthenogenetic species has forced us to 
use parsimony networks, which are well known in the field of 
phylogeograpliy, to demonstrate the peculiarities of the geograph- 
ical distributions of individual clones. Despite the formal similarity 
of the two types of molecular markers, the mechanisms of 
accumulation of nucleotide substitutions (SNPs) in mitochondrial 
DNA molecules are rather different. Here, we used many 
simplifications and assumptions, such as the absence of linkage 
among the three loci, the independence of replication events at 
each of the three loci, the equal rate of microsatellite repeat 
deletions or insertions, and the lack of intergenic recombination 
and genie and chromosomal interference [31]. These processes 
can distort the relationships between the genotypes; hence, in the 
course of parsimony network analysis, we evaluated only 
peculiarities in the geographical distributions of the genotypes. 
Therefore, we were not able to assess the possible time of 
divergence and specify the ancestral clone. 

Low intraspecific variation and a high degree of mtDNA 
similarity regarding the parental species have been a general rule 
in parthenogenetic Darevskia species [6] . These features have been 
explained by young age. Darevsky et al. (198,5) estimated that the 
parthenogenetic Caucasian rock lizards originated after the 
Pleistocene glaciation of the Caucasian mountains, 5000 to 7000 



years ago [16]. Based on the aUozyme analysis of Z). dahli and two 
parental species, Murphy et al. [18] suggested the possibility that 
D. dahli arose in Georgia and dispersed into Armenia. This is 
correlated with today's Georgian localization of the maternal 
species D. mixta and the strong allopatric distribution of parental 
species. It can also mean that the initial microsatellite clones 
detected by us originated in Georgia, and then dispersed to 
Armenia. At the same time, the possibility that these clones 
originated in Armenia could not be excluded, because many 
climatic and vegetative changes have occurred postglacially in the 
Caucasus region, providing an abundance of microhabitats for 
hzard species. An alternative scenario of the origin of D. dahli can 
also be considered. Two populations, Vaagni and Dendropark, 
were characterized by the highest estimates of genetic diversity 
(Table S2). However, in the case of sympatric distribution of the 
parental species in the past, the most probable D. dahli origin 
should be the site of the Vaagni population. The greatest number 
of genotypes, including six specific genotypes (4, 5, and 8-1 1) was 
detected in this population (Figure 3, Table 2). Moreover, the two 
of them, genotypes 10 and 11, originated from independent 
bisexual parental individuals. This scenario supposes the further 
distribution of D. dahli from the Vaagni site into other locations in 
Armenia. Clone 5 might be one of the putative initial clones, 
according to the scheme (Figure 3). Genotypes 4, 2, 9, and 1, 
which are derivatives of clone 5, originated in the course of lizard 
migration in a northern direction (Vaagni to Dendropark), as well 
as in a southern direction (Dzoraget to Papanino, to Fioletovo). In 
the latter case, the river Pambak may be a geographical barrier for 
lizard migrations, which is the reason for the low genetic diversity 
and abundance of clone 1 in the Papanino/Fioletovo populations. 
However, the proposed scenario is not the only one. Regardless of 
where the initial clones were originated, their abundance as well as 
that of clones originated through microsatellite mutations may be 
correlated with a climatic dissimilarity of the habitats [36]. All 
mutant clones should have a later origin relating to the initial 
clones. Among them, under similar ecological conditions, more 
abundant clones probably arose first. 

The results obtained for population genetic parameters may be 
useful for comparing the samples studied with the parental species. 
Moreover, the estimates of population genetic indices may be 
useful when comparing parthenogenetic reptiles belonging to 
different genera. 

It should be noted that our microsatellite genotyping data, as 
well as allozyme [18] and mtDNA [8] stucUes, failed to cUstinguish 
the two color varieties ("yellow" and "common" morphological 
variants) of D. dahli. Among the 1 1 1 indi\'iduals of Z). dahli, only 14 
belong to the "yellow" morphological variant. All these lizards 
were found in the Papanino population and had genotype 1. We 
believe that two morphological variants of D. dahli are genetically 
determined because they also differ functionally, namely in egg 
production. Further research using new microsatellite loci will 
probably differentiate these forms. 

In summary, our data suggest that the clonal divc-rsity observed 
for Z). dahli derives from both multiple interspecific hybridization 
origins and postformation microsatellite mutations of the initial 
clones. The methodological approach developed in this study, 
which was based on the detection of genotype-specific markers 
formed by combinations of parental-specific markers, can be 
further applied for the elucidation of the origin and evolution of 
unisexual species. 
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Supporting Information 

Table SI Polymorphisms of three microsatellite loci 
Du215, Du281, and Du323 for nine genotypes in /). dahli 
(see Figure 2). GATA repeat is designated by the letter A. 

Numbers before the letter A designate the number of the repeat on 
each of the two chromosomes that are separated in the Table by 
forward slash "/". In concatenated sequences used for parsimonies 
network reconstruction, the absence of a GATA repeat is 
designated by dash (— ). 
(DOC) 
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